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Abstract 

We present a method for computing the PDF and CDF of a non-negative infinitely divisible random 
variable X. Our method uses the Levy-Khintchine representation of the Laplace transform Ee~ — 
e~'^^'^\ where (j) is the Laplace exponent. We apply the Post-Widder method for Laplace transform 
^^ ' inversion combined with a sequence convergence accelerator to obtain accurate results. We demonstrate 

^ , , this technique on several examples including the stable distribution, mixtures thereof, and integrals with 

respect to non-negative Levy processes. Software written to implement this method is available from the 
authors and we illustrate its use at the end of the paper. 



1 Introduction 



Let X be a non-negative random variable. The distribution of X is said to be infinitely divisible (ID) if for 
^ I any positive integer n, we can find i.i.d. random variables X^ „, i = 1,2, . . . ,n such that 



X — Xi^n + X2 



^^ , For background, see |22] . [3]. There are many examples of such distributions, including the gamma distribu- 

10 ' tion, compound Poisson distributions, inverse Gaussian distribution and right-skewed stable distributions. 

^^ , These distributions are also central in the study of non-decreasing Levy processes. In this paper, we give 

a method for numerically computing the probability density function (PDF) and cumulative distribution 

function (CDF) of a non-negative infinitely divisible random variable. 

Our starting point is the Lcvy-Khintchine (LK) formula ([i], IS],[1T]), which in the case of non-negative 

ID random variables states that the Laplace transform of X, 

^'- i^(A) = Ee-^^ = e-^(^), A > 0, (1) 

has an exponent (t){X) which is called the Laplace exponent, which can be written as 

0(A) =aX+ (1 - e-^")n(du). (2) 

Jo 

Here, a > is a shift and 11 is a measure on (0, 00) which satisfies 

(1 Aa;)n(da;) <oo. (3) 
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Since we are interested primarily in PDFs and CDFs, wc will assume that a — throughout, since a positive 
a only shifts the PDF/CDF. Let fx{x) and Fx{x) denote the PDF and CDF of X, respectively (we will 
promptly drop the subscript X on these functions when it is clear from the context). The Laplace transforms 
of the PDF and CDF can be easily obtained from the LK formula (we will use a tilde to denote a Laplace 
transform) : 

V'(A) = /(A) = / e-^^f{x)dx, 
Jo 

and a simple application of Fubini's theorem implies 

^{X) ^ F{\) = ^°° e-^- (^£ f{y)dy^ '^''^kj^ ^-'''f{y)dv = ^. (4) 

Thus, obtaining / and i^ is a matter of inverting a Laplace transform. Generally, this task is not easy. 
Typically it is done by complex integration of the Laplace transform (see [2]), which can be difficult if the 
integrands are slowly-decaying, oscillatory functions. This causes many numerical integration methods to 
converge slowly. Here, we apply a different method of Laplace inversion known as the Post-Widder (PW) 
method ([T], Theorem 2 or [7], section VIL6). It is based on the fact that under weak conditions on a function 
/, we have 

where /^'^"^^(fc/x) denotes the (fc — 1)*'' derivative of / evaluated at k/x. Thus, instead integrating the 
Laplace transform, wc are taking arbitrarily high derivatives. 

The obvious challenge in using this method is computing high derivatives of / in ([S]) to approximate the 
limit. Some methods have been developed to do this in general (see for instance, [10] or [2]), however most 
involve complex integration. The method we describe here for the case of infinitely divisible distributions is 
both easy to implement, and only involves (at worst) integration of a real valued non-negative exponentially 
decaying function. The method we use here combines both computation of a finite number of terms of the 
sequence in ([5|) and an numerical extrapolation method to approximate the limit, [5]. A similar method was 
used to compute mean first-passage times of Levy subordinators in |23j . 

This paper is organized as follows. In Section [2] we give an overview of the approximation method used 
based on the (PW) formula. In Section [3l we outline the algorithm for finding / and F, as well as numerical 
issues that may arise in the computation. We test the method in cases where the PDF and CDF are known 
in closed form in Section U) In Section [S] we apply our method to a collection of examples. The software 
for implementing the methods described here is freely available form the authors and its use it described in 
Section [6] 

2 Post-Widder method with extrapolation 

Given a continuous function g{x), a; > 0, which is bounded as a; — > oo, we will denote its fc*'' PW approxi- 
mation as 

(—1) I k\ ~(k-l) (^ 



The fact that gk(x) — )■ g{x) as fc — > oo can be seen by approximating g(x) by £(7(1^.), where Yj, = k^^{Yi + 
i^2 + ' • • + i^fe), and the l^'s are i.i.d. gamma random variables with mean x and variance x^ and then applying 
the law of large numbers (see Section VII. 6 in [7j). 

The convergence of gk{x) to g{x) is slow in general. To illustrate this point, let g be a PDF / of an 
inverse Gaussian distribution for which the Laplace transform is given by /(A) = cxp(— A^^") (in this case 



the density f{x) is known in closed form, see Section 0]). In Figure [U we plot the exact formula for / as well 
as /fe for A: = 1, 10 & 50. Notice that even with 49 derivatives of /, the approximation is still poor. In fact, 
this is a general feature of the PW formula, as it has been shown ([10]) that the errors efc(/; x) = f{x) — fk{x) 
have a power series expansion: 



ikif;x) 



m— 1 



(7) 



where the coefficients amiin) are given by 



am{x) =^f 



(j+m) 



{x)x- 



j+m 



J = l 



d(j + m, m) 
(j+m)! ' 



where d{j + m,m) arc the associated Stirling numbers of the first kind ([T^, Chapter 4, section 4). Thus, 
the convergence of fk is in general 0{l/k) corresponding to the m = 1 term. This method alone is thus 
inadequate for computing the inverse Laplace transform to a high level of precision. 




Figure 1: The exact form of / together with the PW approximations /i, /lo and /50 

To obtain high precision, it is necessary to couple the PW approximations with a convergence acceleration 
method which extrapolates the limit in ([5|) based on a finite collection of terms in the sequence. Thus, our 
method for computing the PDF / at a; > involves two components: 

1. For some sequence ki < k2 < ■ ■ ■ < k^ , compute the approximates fki{x). 

2. Use the points {ki,fk^{x)), . . . , ikN,fk^{x)) to extrapolate limfe^oo fkix). 

And for the CDF, replace / with F above. We focus first on step 1, which in sight of ^ involves computing 
high derivatives of ip or ^. Next, we'll discuss two methods of convergence acceleration to address step 2. 

Remark: Our method produces the best results if the PDF and/or CDF are smooth functions for a; > 0. 
That being said, this method still produces useful results in the non-smooth case, however it may fail to 



converge near points where the function lacks smoothness. This method should not be used for distributions 
which contain atoms, such as the Poisson distribution, or even distributions whose density lacks smoothness, 
such as a compound Poisson distributions with bounded jump distribution. For such distributions another 
method based on numerically solving a Kolmogorov- Feller forward equation, is more suitable, see |24j . 



2.1 Computing derivatives 

Derivatives of ^ 

We begin by presenting two methods of computing derivatives of ?/'(A) = exp(— (/)(A)). The first method is 
based on the simple observation that 

^'(A) = -f (A)V(A). 

Thus, Leibnitz's formula implies that for any A: > 1, the {k + 1)*'' derivative of ^ is given by 

k 






(8) 



Notice that for A fixed, the values of tp{X), V''(A), . . . , V''-'^-' can be computed recursively using (|S]) if one has 
first computed 4>^^^ (A) for j = 1, 2, . . . , fc + 1. 

Alternatively, the fc*'' derivative of tp{X) = exp{—(f>{X)) can be given directly in terms of the derivatives 
of (j) using Faa di Bruno's formula ([18|, Chapter 2, section 8), which yields in our case 



^W(A) 



dX^' 



-./.(A) _ ~4>W 



E 



fc! 



mi\m2\ ■ ■ ■ ruk^. 



>'^HA) 
1! 



-0W(A)Y 

ki ; 



where the sum runs over all ttii, . . . ,mk such that mi + 2m2 + • • • + knik = k. This can expressed more 
simply in terms of the /c*'' complete Bell polynomial, Bk{xi, X2, ■ ■ ■ ,Xk), as 



^W (A) = e-^(^)i3fe(-(/)(i) (A), -0(2) (A), . . . , .^W (a)). 
See [17] or [19]. Furthermore, Bk can be given in terms of the following k x k determinant: 



(9) 



Bk{xi,...,Xk) 



XI {\')X2 C=-^)X3 {'-')X, 

-1 XI {\'')X2 {^f)x, 

-1 xi '^-^^ 

0-1 



Xl 

















Xk 




Xk-l 




Xk-2 




Xk-3 


-1 


Xl 



(10) 



This provides another way for computing ip^'^' . 

Both methods above require the derivatives of the Laplace exponent (f>. Fortunately, these can be com- 
puted directly, since the LK formula implies that for any A > 0, 



(1 - e--'")n(ffu). 



■')(A) = <; 



(-1) 



n+l 



u"e-^" 



U{du), 



n = 



n> 1 



(11) 



Notice that passing derivatives through the integrals above is justified by the integrabihty assumption ^, 
since for any n > 1, x"e~^^ < C„.a(1 A x) uniformly on {x > 0} for some constant C„.a > 0. For 
many examples, these integrals have a closed form expression and can be computed easily. If a closed form 
is not available, many numerical integration methods are effective as these integrands are non-negative, 
exponentially decaying functions. In this case, it is important to use a small relative error tolerance, as 
(/)*^")(A) can become extremely small for large A. 

Derivatives of ^I* 

For \I' defined in Q, we again apply Leibnitz's formula to ^ and obtain 



$W(A)=f('=)(A) 



d'' Vj(A) 
dX'' A 



E 



^^i-iy^i^^'-'H^i 



X3 + 1 



When combined with the Post-Widder formula ([6]) , we see a good deal of cancellation: 



(fc-1)! \x 



(-1) 



fc-1 



(fc-1)! 



Ef':')(-i)^ 



(_i)fc+j-i 



fc-1 



3=0 
X 



j! 



J y ■ ■ (k/xy 



-^ 



(fc-i-j) 



(k/x) 



(12) 



Note that this requires the values ^{X),ip'{X), . . . , V'^"-' (A), which can be computed using (|S]) and ([TT|) hence 
the recursive formula ([5]) is more convenient here than ^ since it computes all derivatives of ip. 



Other useful derivatives 

In some applications, the derivatives of the PDF are also of interest. Assuming the first g — 1 derivatives of 
/ vanish at 0, i.e. /(O) = /'(O) = • • • = /(«-i)(0) = 0, then /^''''{x) has Laplace transform 

/(^(A) = AV(A) - A«- V(0) - A'- V'(0) r~\0) 

= AV(A), 

where tjj{X) — /(A). In which case, Leibnitz's formula gives 

gA(fc-l) 



Jfc ^fc 



J=0 



«! 



i/ (9-j)! 



A'3-J^('=-J)(A) 



(13) 



Thus, (HI) and (|lip can be applied to compute the derivatives of /. 



2.2 Extrapolation 

Wc want to compute the PDF f{x) or CDF F{x) by approximating the limit as fc —> od in the Post Widder 
formula ([S]). By letting h = fc^^, ([7]) implies that f{x) can be written as 

oo 

fix) = Mx) + ^ a,„/i™ (14) 

771—1 

We shall truncate the error series X^tti^i at m = A^ — 1 and consider non-consecutive fci < fc2 < ■ • • < fcn- If 
hi — k^ , i ~ 1, . . . , N, we can use the points {hi, fk^ (x)), (/12, fk2 (2;)), ■ • • , {hN, A-n (2^)) to estimate the first 
A^ — 1 unknown coefBcients a„i. Indeed, the (approximate ) N x N system 

fkAx)~I{x)-aih,-a2h] anh^-\ .? = 1, 2, . . . ,iV (15) 

can be solved for the N unknowns f{x),ai, . . . ,a7v-i. Since we are only interested in f{x), it is enough to 
only compute the first row of the inverse of the matrix corresponding to the system (|15p . This is essentially 
done in the polynomial extrapolation method described below in order to obtain an accurate approximation 
of the limit Imik^oo fk{x). There is a vast literature on such extrapolation methods and the errors associated 
to them, see [TT] for a review. 

Below, we briefly review two techniques for approximating ,f{x) also used in [8]. We call these methods 
of extrapolation. The first is based on polynomial interpolation, and the second based on rational function 
(Pade) interpolation. Not surprisingly, the rational extrapolation provides faster convergence in many cases, 
however wc found it more susceptible to numerical instability for larger N and fc^r. Therefore, the "better" 
choice depends on the particular example and the desired accuracy. 

Polynomial extrapolation 

Given the points {hi,fk^{x)), (/12, ^2(2;)), . . • , {h^, fk^ix)), the iV - 1 degree Lagrange polynomial PNix) 
which passes through these points is given by 

1=1 \j^i * 



Thus, taking h = Q (or, k = 00) in this polynomial provides the approximation 

- h, _(-l)^-ifcf-i 

h 



m^pm^Y..hM. wh™ ,.n^ = t^^. 



Wc thus use a linear combination of the approximations fk{x) to obtain a more accurate approximation of 

Rational extrapolation 

As an alternative to polynomial extrapolation, on may instead fit a rational function to the points (/ii, fki{x)), 
. . . ,{hN, /Arfc7v(a;)) of the form 



RnW = 



P,ih) 



QAhY 

where P^{h) = J2i=iPi^^ i^ ^ polynomial of degree /i ~ [N/2\ and Qv{h) = Yli^QQih'^ is a polynomial of 
degree v ^ N — \_N/2\ . Here, the coefficients p^, qi are chosen so RN{hi) ~ fki{x)- When implementing this 



method, these coefheients are not computed directly hke in the polynomial extrapolation case, but instead 
Rn{0) is computed iteratively in a triangular array using the following recursive formulas (see [11], section 
13) 

BL = Kt\ + -r/ "-\ where Ql ^—^ ^—^ 



i^i^m l^i+m ^m-1 ^ ^)ri- 



2 



Wc then take /(x) w i?]y = Rn{0) as our approximation since /i = 1/k -^ corresponds to fc ^- cx). Notice 
that unlike the polynomial extrapolation, this transformation is nonlinear in the fki (x) . 

Error bounds 

One major drawback to using extrapolation techniques is that the error made in the approximation is difficult 
to bound analytically (see, for instance. Theorem 3 in [6] and equation (13) in [8]). One can however, obtain 
an asymptotic error bound for a rational or polynomial approximation. Let hi < h2 < hs . . . be a sequence 
for which the approximant fki{x) can be computed with ki = l/hi. Also, let Pjv(O) be the approximation of 
/(x) obtained by using polynomial extrapolation with /ii, /12, . . . , /ijv, with N > 2 (all the following formulas 
also hold with rational extrapolation by replacing P/v(0) with i?7v(0)). In |Q, Bulirsch and Stoer (BS) 
construct a second estimate of /(x), PAr(O), with the property that 

lim £^(t/M ^ _,. (16) 

^^-Pw(0)-/(x) 

Thus, asymptotically, P/v(0) is as good an approximation to f{x) as PAr(O), except that it approaches f{x) 
from the opposite direction. This second approximation is obtained using a linear combination of the form 

Pjv(0)-(l + a)PAr+i(0)-aPw(0), (17) 

and choosing the constant a to change the sign of the leading term in the error Pn{0) — f{x)- BS show this 
constant a is given by 

/lAf + l 

This allows us to construct a numerical bound on the relative error made in this method. From (|16p . we 
have _ 

Pn{0)-Pn{0) ^ Pn{Q) -J{x) ^1^2 (19) 

/(x)-PAr(O) /(x)-Pjv(O) 

Fix any -q S (5,1)- From p^ and since ?7 > ^, there exists NQ{r]) large enough such that for N > Nq^tj), 
\f{x) — P/v(0)| < 77|PAr(0) — PAr(0)|. Using a similar argument, there also exists Mo{r]) large enough such 
that for N > Moiv), l-FW(O) - /(x)| < ?7|PAr(0) - Pv(0)|. Moreover, since PAr(O) -)■ f{x) and r? < 1, we can 
pick Ni{ri) large enough such that |Pjv(0)//(x)| < 1/?? for aU N > Ni{ri). 

We thus can refine our estimate as f{x) « (Pjv(O) + Pv(0))/2. If A^ > max(iVo(f?), ^0(77)), 



PNiO)+PN{0) _ ^(^) 



^\jm^m^\m^iEm<,\P,^o)^P.io)i<iP.io)^PM 



and if iV > mi,xiNo{v),Moiv),Ni{r^)), 



Pw(0)+Piv(0) 
2 



/(^) 



/(^) 



< 



< 



< 



\PN{0)-.f{x)\ |/(x)-Pa.(0)| 
2f{x) ^ 2/(x) 

Pn{0)\ \Pn{0)-Pn{0)\ 



fix) J Pn{0) 

|Pjv(0)--Pw(0)| 
Pn{0) 



Therefore, 



Theorem 2.1 Suppose f{x) ^ and let PAr(O) &e an approximation of f{x) obtained using either polynomial 
extrapolation or rational extrapolation, and let P/v(0) be defined as in |i7[). For N large enough, we have 



Pn{0) + Pn{0) 



-f{^) 



< |PAr(0)-Pjv(0)| ano 



Pw(0)+Pw(0) 



/(^) 



/(-^) 



< 



|Pjv(0)-PAr(0)| 

Piv(O) 



(20) 



where the right hand sides converge to as N ^- oo. 



This provides a numerical bound on the absolute and relative errors of the approximation and suggests 
increasing TV until either \Pn{0) — Pn{0)\ or |P;v(0) — Pn{0)\/Pn{0) is smaller than a prescribed value. 

3 Implementation 

We begin with a description of the algorithm for computing the PDF and CDF based using the recursive 
formula (|8]) and the polynomial extrapolation method. This method is implemented in MATLAg^ and 
Mathematicq^ The method using rational extrapolation is also an option in the MATLAB implementation, 
but is not described here. 



Algorithm for computing the PDF / at a; > 



1. Choose the sequence 1 = fco < fci < ^2 < • ■ ■ < fcw and a relative error tolerance e > 0. We found 
kj = lOj, with N lying around 8 to be effective for e = 10~^. Initialize the following arrays: 

Variable Size Purpose 

D<l) (fcjv — 1) X A^ Holds derivatives of </> 

D.0 k^ X N Holds derivatives of ip 

P I X N Holds Post-Widder approximations of / 

/ 1 X A^ Holds extrapolated approximations of / 

Set j = 1. 

2. Compute {D.^)ij as 



^MATLAB is computational software package developed by The Mathworks 
■^Mathematica is a computational software package developed by Wolfram Research 



3. For i — fcj-i, . . . , kj — 1 , compute the r row of D^ as 

|0 otherwise ' 

(^''^ is computed using the integrals in (|lip . The (i + 1)*'' row of D^, is then computed as 

[ otherwise 

4. Compute the j*'' Post-Widder approximation: 

^^" - TfcT^T)! UJ ^ ''^'^''' 

5. If j = 1, set /i = Pi, J = 2, and go to step 3. Otherwise, compute the j*'* extrapolation as 

fj = y^CrPr, where Cr = TT 






Set a = 1 + 2(fcj/fci + l)-i. If 1(1 + a)(/j - fj-i)\/fj < e or if j = TV, return /j. Otherwise, set 
j = j + 1 and go to step 3. 

The method for computing the CDF is similar, except that ([T2|) is used instead for step 4. 



Algorithm for computing the CDF F at x > 



Follow steps for PDF computation above except replace step 4. with 
4.* Compute the j*'* Post-Widder approximation: 

^^- - ^^ (fc-373T)I [-) (^'>)^-- 

We'll conclude this section with a collection of remarks regarding implementing this procedure. 

Remarks 

1. Binomial coefficients are reused multiple times in step 3. We found it useful to compute and store 
rows of Pascal's triangle as needed. 

2. In MATLAB, all computations above can be "vectorized" to maximize speed. 

3. It is essential to compute the derivatives of (fi to as many significant digits of accuracy in step 3 as 
you want in your final result. In many examples, a closed form for (/>'■"■' can be found. See Section 
[5] for examples where these integrals are computed without a closed form. 



4. To avoid overflow, ratios of large numbers like those seen in steps 4. and 4.* should be computed 
to incorporate reduction. For example, using identities such as 



(fc, - 1)! 



exp I kjlog{kj/x) - Y^ log(j) 



will alleviate overflow. 

In double precision arithmetic, it is often impossible to take fc^r much bigger than 100 as un- 
derflow, overflow and numerical instability become unavoidable. In most cases, adequate con- 
vergence is met much before such large k values are needed, however if one must go further, 
it will become necessary to work with a high-precision arithmetic package. Mathematical for 
example, has powerful multiprecision capabilities. If one requires more speed, a well devel- 
oped and documented library for multiple precision arithmetic for C++/Fortran is available at 
http: //crd. Ibl .govZ-dhbailey/mpdist ■[ 

When the true value of / is very close to (/(x) <^ 10^®) the extrapolation procedure might 
"overshoot" and return a negative value for the density. In this case, is a good approximation 
as you can be sure that the PDF or CDF takes an extremely small value in this case. And similarly 
when the CDF takes a value greater than 1. 



4 Testing the method 

In this section we use Mathematica and consider two examples for which it is possible to compute the PDF 
and CDF exactly. The first example is the chi-squared distribution x^ with one degree of freedom. For this 
distribution, 

/,.(x)^^.-V2,-./2^ F,.(.)=erf(^y|) 
^^.(A)=.4.(A) = (l + 2A)-i/2 



0^2 (A) = -log ^^2 (A) = -log(l + 2A) 



(1 



-Xu\ 



-«/2 



2u 



du. 



'^5^(A) = (-1) 



n-t-l 



n —An 
U e 



,-«/2 



2u 



du 



'-""'f "„•-.-»■/'».„.'-""■"- "74 + A 



n > 1. 



Results over a range of inputs is shown in Table [TJ We show the a;- values considered as well as the value 
of N required to obtain a relative precision of 10^^ and 10^^^. In Figure [21 we also plot the relative error 
as a function of N . Notice that the value of N required is substantially higher when the value of the PDF 
or CDF takes on very small values, in particular, when /;^2(50) « 7.83 x 10~^^. This is because /fc(x) over 
approximates / at these points. 



10 




Relative Error 



5 10 15 

(a) Relative error for PDF of x^ distribution 




5 10 15 

(b) Relative error for CDF of x^ distribution 



Figure 2: Plot of relative error of our method versus the truncation level N for the PDF and CDF of the 
X^ distribution at three values of x. 




5 10 15 20 

(a) Relative error for PDF of IG distribution. 



Relative Error 




5 10 15 20 

(b) Relative error for CDF of IG distribution 



Figure 3: Plot of relative error of our method verses the truncation level N for the PDF and CDF of the 
IG distribution at three values of x. 

The second example we considered is the inverse Gaussian (IG) distribution (p], Example 1.3.21), for 
which 



/. 



IG 



(x) = ^=L=e-i/(4-), FjGix) = l-cri(^) 
Vinx^ \'s/2xj 



^igW = fiaW = e 



-V\ 



biaiX) = - logV^/G(A) = VA = / (1 - e-^") 

Jo 



oo 



U 



-3/2 



2V^ 



du 



<^fe'(A)=^"'^ 



71+1 /-OO 



2V^ 



n-3/2 -Am J ( 1)" j^/ i/o\\l/2-n ( 1)" (^^ •^)-^l/2-n 

J. ' e du= — ——j^—i (n — 1/2)A ' = ——- — -r- rr;— A ' 

20F ^ ' ' 22("-i)(n-2)! 



Table [T] also shows the values of x considered, the value of ficix), and the value of N needed to obtain 
relative errors of 10^^ and 10^^^. Figure [3] shows how the relative error behaves as a function of A'' for three 
values of x. Similarly to the x^ distribution, N is largest when the value of the PDF is very close to 0. 



11 



Chi-squared 


X 


10-^ 


10-^ 


1 


10 


20 


50 


f^2 {x) 


1.26 X 10^ 


1.20 


0.242 


8.50 X 10-* 


4.05 X 10-"^ 


7.83 X 10-^^ 


Ne 


4 


4 


4 


6 


9 


14 


Ni5 


8 


8 


9 


12 


15 


21 


FA^) 


2.52 X 10-^ 


0.248 


0.683 


0.998 


1-7.7X 10-« 


1 -1.54X 10-" 


m 


3 


3 


3 


4 


4 


1 


Ni5 


8 


8 


8 


11 


11 


12 


Inverse Gaussian 


X 


0.01 


0.02 


0.1 


1 


100 


1000 


ficix) 


3.92 X 10-« 


3.71 X 10-* 


0.732 


0.220 


2.81 X 10-" 


8.92 X 10-^ 


Ne 


12 


8 


6 


5 


6 


6 


Ni5 


21 


17 


14 


11 


10 


11 


Fig{x) 


1.53 X 10"^^ 


5.733 X 10-' 


0.025 


0.480 


0.944 


0.982 


Ne 


12 


9 


6 


4 


4 


3 


Ni5 


21 


17 


13 


10 


8 


8 



Table 1: Using ki — lOi, this table shows the value of N required to obtain 6 digits of accuracy (Nq) and 
15 digits of accuracy (A^is) at the given values of x using polynomial interpolation. Notice that it becomes 
harder to approximate these functions to a small relative error when they become very small. 

5 Applications to ID distributions which are not known in closed 
form 



Wc will now consider examples of non-negative ID distributions for which the PDF and CDF are not known 
in closed form. We will compute them numerically using our method. The resulting plots of the PDFs 
and CDFs are obtained using MATLAB. Each plot was generated in about 1 second. In order to apply 
our method, we much first write the Laplace transform in LK form, and then find 0(A) and the derivatives 
(?!)*^"^ (A) by computing the integrals in ([TT|). In almost every case below, these integrals can be given in "closed 
form" , which is to say they can at least be written in terms of special functions for which there are efficient 
methods for computation. We will assume through out that the drift a = in ^. 

The following special functions will appear throughout this section, we provide their definitions here for 
convenience: 



Gamma function 



Lower incomplete gamma function 



Upper incomplete gamma function 



Entire exponential integral 



Dilogatithm 



r(a)=/ z°-^e-^dz, a ^ 0, -1, -2, . . . 
Jo 



T{a,b)= z'^^'e-^'dz, 6 > 0, a G 

Jb 

r I- er"" 

Ein(a) = / dz, a S R 

Jo z 

L,ia)^ r^-^^dz, a>0 
Ji z-1 



See 



chapters 25, 37, 43 and 45 for more discussion of these functions and methods for efficient compu- 
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"-1 5" 



tation. For an integer ri > 1, we have T{n) — [n — 1)! and 

7(n,6) = (n-l)!(l-e-''^^), (21) 

V m=0 ^^' / 

see [5], Equation 3.351.1. 

5.1 Right-skewed stable distributions 

Here we consider the class of examples with Laplace exponent given by 



0(A) = / \Pp{dl3) = / (1 - e-^") 
Jo Jo 



1 u-^-^ 



p{d(i) 



du, (22) 



/o r(-/3) 

where p is a measure supported on (0, 1). The Levy measure 11 has a density is given by 

Mixtures of this form are considered in [13] and [12], in which these distributions are used in models of 
anomalous diffusion. 

The Laplace exponent is expressed in p2p and the derivatives 0^"-* for n > 1 can be computed by using 
(jlip . a change in the order of integration, and the definition of the gamma function: 



•>(A) ^ (-!,-« f .-»- (/' ^p,.«) .„ . (-1,- / ^q^pm 



n> 1. 



n-l 



Since n takes only integer values, r(n — /3) ~ r(— /?) nm=o('^~/^)' ^^"^ ^^ ^^^ above can be simplified further 



as 



(„)(^) = (-1)"+^ /' ^!ri(!i^p(d/3) = (-1)"+^ /' A^-" n (™ - /3)Krf/3) 

^0 n-P) Jo „^o 

1 " r^ 

= -^^E^'-W^n'"^ (23) 



m— 1 



where Cm(A) = L [3"^ X^ p{d[3) and {S*,!™ }, ri>0, 0<m<n are the Stirling numbers of the first kind ("IIS]. 

page 162). These are such that riJ^^o(/3 — m) = X)m=i SnP"^, and can computed with a triangular array 
similarly to Pascal's triangle using the recursion formula 

5^ = 1, 5(°)=0, n>l 
St^=st-,'^-{n-l)stl r.,m>l. 
Let us now consider special cases of p (below, 5 denotes the dirac (5-distribution). 
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5.1.1 Right-skewed a-stable distributions: p{df3) — S{f3 — a)d/3, < a < 1 

Indeed, an important example of this distribution is tlie right skewed a-stable distributions, for wliiclr p is a 
point mass at /3 = a, with < a < 1 (using Proposition 1.2.11 in (5D], this distribution corresponds to the 
stable distribution S{a, 1, cos(7ra/2)^/", 0)). These distributions lie in the family of scaling limits for sums 
of non-negative i.i.d. random variables with infinite mean. Form (j23p . 0^"^ can be computed in various ways: 

U^j = (-1) Ti-a) "^ J_J_(a-m) = -A 2^ a Si, > 

The PDFs and CDFs of the right skewed a-stable distribution are plotted using our method in Figure |4] for 
several values of a. For example, we compute the PDF of a (l/2)-stable distribution at x = 1 and obtain 
(in Mathematica) 

/(I) = 0.219695644733861 

This is exact to 15 decimal places and the computation took approximately half a second. An alternative 
method for this case is given in [T4] . 



5.1.2 Sums of right-skewed a-stable distributions: p{d(3) = 2^djS{(3~aj)d(3, dj > 0, /^ rf j = 1 

The previous case can easily be generalized to weighted sums of independent a-stable random variables. 
From (1^51) it follows that in this case. 



j = l ^ ^' j = l m=0 j=l \m=l / 

Plot of PDFs of such distributions can be see in Figure [SJ 

5.1.3 A "uniform mixture" of a-stable distributions: p{d/3) = dfi, j3 e (0,1) 

This is an example of a distribution with no finite moments. The Laplace exponent (j22p is given by 

0(A) = /' X^dp = t e/i°g^d/3 = ^^^ if A ^ 1, 
Jo Jo log A 

and is 1 if A = 1. Higher derivatives of (f> can be computed using (|23p. The coefficients c„i(A), m — 1,2, . . . ,n 
can be computed in a few different ways. First, if log A < 0, we have 

c™(A) = [' P-X^dp = [' re^^^^^dp = "^;V"\;J:f^ (24) 

where 7(a, 6) is the lower incomplete gamma function (see above). The formula \2A\ can still be used for 
log A > 0, however this requires analytic continuation of 7 which is not always easy to compute. 

As another approach, we compute Cm(A) by treating the cases | log A| < 1 and > 1 separately. If | log A| < 1 
notice that 

(logA)J 



c™(A)=/ /3'"ei°s(^)^d/3 = ^(logA)W ^d/? ^ ^ 

"'0 ,-n "'o ■?■ o-n 



— ' {m -|- 7 -t- 1)7! 
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(25) 



If I logAj < 1; taking the first 18 terms in this series gives an absolute error of < 10 ^^. For | log A| > 1, the 
c,„'s can be computed recursively by applying integration by parts: 



CO (A) 

Cm (A) 



,log(A)0^^ ^ 



A-1 

bi(A) 



log(A) 



Since this procedure involves a division by log A, it should only be used when | log A| > 1. 



a-stable PDFs 



Q-stable CDFs 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 




1.2 1.4 1.6 



Figure 4: a-stable PDFs (left) and CDFs (right) for a = i/10 for i = 1, 2, . . . , 9. The PDF for a = 0.9 
corresponds to the right-most peak, and a = .1 corresponds to the left-most peak (which isn't visible). The 
CDF with the steepest slope around x = 0.8 corresponds to a = 0.9. 



5.2 Integrals of non-random functions with respect to a Poisson random mea- 
sure 

In this section we consider the large class of non-negative ID distributions which can be expressed in terms 
of a Poisson stochastic integral of a non random kernel. For a recent review of these integrals, see |16j . 

Let (r^,^) be a probability space let B denote the Borel sigma field on R". Let iV(-) denote a indepen- 
dently scattered Poisson random measure with control measure /z, that is, ^ is a measure on B and TV is a 
function 17 x 6 ^ Z+ U {0} such that 

(i) N{A) _L N{B) if A, B G S are disjoint. 

n n 

(ii) N{ \^ Ai)^^ N{Ai) if A, are disjoint. 



m— 1 



?n — 1 



15 



1.2 



0.6 



0.2 




m = 1 
m = 2 

TO = 3 

TO, = 4 
TO, = 5 

Uniform Mix 



0.25 0.5 0.75 1 1.25 1.5 1.75 2 




1.25 1.5 



Figure 5: (Left figure) This is the plot of the PDF corresponding to the choice Puj{du) = (1 — uj)S{x — 
0.4) + uj6{x — 0.8) for various choices of < <^ < 1- By increasing w, the PDF for a = 0.4 (on the left) is 
morphing into the PDF for a = 0.8 (on the right). (Right figure) These are PDFs corresponding to the sums 

Prn{du) = — X^I^Li '^(^ Tt) ^'''^ "^ ~ 1, 2, . . . , 5 together with the uniform continuous mixture on (0, 1), 

which has the highest peak. 



(iii) For each A E B, N{A) has a Poisson distribution with mean ii{A). 
Given such a pair {N, /i) one can define the following stochastic integral 



Ha) 



g{x)N{dx) 



(26) 



for a function g on R", for which J^^ nim{l,g{x))^{dx) < oo and, for our purposes, is non-negative. In this 
case, the random variable I{g) is also non-ncgativc and has Laplace transform 






=.-As(s) 



)Kds) 



A>0. 



(27) 



In order to compute the PDF and CDF of /(g) using our method, (P7|) must first be rewritten in 
LK form. In many cases, this can be done with some suitable change of variables u ~ ^(s) satisfying 
ui = (4>(s))i = g{s). In this case, (P?]) becomes 



Ee-^^(») = exp (- f (1 - e^^"i)ng(dui)' 
where for any Borel set A C M+ , the Levy measure Ilg {A) can be expressed formally as 



^M) 



AnM"-i 



|J(ui,U2,...,w„)|(mo $ ^){dui,du2,...,dun), 



(28) 



(29) 
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where \J\ is the Jacobian d^~^(ui, . . . , Un)/d{ui, . . . , u„). 

To illustrate this, let us now focus on special cases, in dimensions n ~1 and n ~2, where this change of 
variables can be made and our method applied. 

5.2.1 One dimensional Poisson integral 

Assume n — \ and that the integrand g is a monotone, no n- negative function with inverse g^^ ■ In this case, 
(j27p can be rewritten in LK form using the change of variables u — g{s): 



Ee--^^(9) = exp ( - / (1 - e--^")ng(dM) 
Since u = g'^^{g{u)), we get 1 = {g' o g~^){u) (g^^)'(w) and hence, 

n.(rf«) = ^f°--))(") (^og-i)(d.)), (30) 

lig'og l)(u)| 
where g{{0, oo)) denotes the image of (0, oo) under g. Suppose that we want to get the PDF and CDF of 

/•OO 

lig) = / e-'^mids) 



where g{s) = e *''', r; > is a parameter and the control measure /i is Lebesgue. In this case, a simple 
calculation shows that (l30l) becomes 



Ugidu) = \^-j 1(0,1] (u)^^. 

Thus, for this example, 

0(A) = / (1 - e"^") (-) du = ?7Ein(A), 

where Ein(A) is the entire exponential integral defined earlier. The derivatives 0'^")(A), for n > 1 can be 
given in closed form 

0(")(A) = (-1)"+^ fu-'e-^-du = I^H;;^^(n,A) = ^^^ f(n- 1)! - e"^ g ^^^"l ' 

where the last equality follows from ([21]). We've plotted the PDF and CDF of the random variable /(e^'*/'') 
for various values of 77 in Figure ([6]) using our method. Note that since the range of integration here is (0, cxd), 
the method described in [21] doesn't readily apply. 

5.2.2 Integration with respect to non-negative Levy process with a non-negative kernel 

In this example, we generalize the previous case by looking at integration with respect to a non-negative 
Levy process, or equivalently, a one-dimensional non-negative ID random measure L with control measure 
^. L is a random measure which satisfies the same conditions as the Poisson random measure iV, except 
condition (iii) is replaced by 
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1.5 
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0.75 
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0.25 









7? - 0.5 

r? = 0.9 

r? = l _ 

7? = 1.1 

77 = 2 


i^; 





•,/^Qs^ 


^..^^ 



1 

0.9 
0.8 
0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 




il^ 


' 


\[/ 


r? = 0.5 
?7 - 0.9 

7?=1 

r? = 1.1 - 

7? = 2 
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1.5 



2.5 



3.5 



Figure 6: Plots of the PDF and CDF for the Poisson stochastic integral /(e "/'') with various values of the 
parameter 77, which plays the role of a shape parameter. The PDF for r\ = 0.5 is the highest on the left. 



(iii)' There exists a Levy measure IT such that for any A G B, the distribution of L{A) has Laplace transform 

Eg-AL(A) ^ g^p (~fi{A) / (1 - e-^'')n{du)) . 

Notice that the Poisson random measure N corresponds to the choice Il{du) = 6{u — l)du. For a given 
function g, the stochastic integral I hid) can be defined in terms of a two-dimensional Poisson stochastic 
integral: 

/•OO /"OO /"OO 

Il{9)= g{s)L{ds)= / ug{z)N{du,dz), (31) 

Jo Jo Jo 

where the control measure of N is given now by Ii{du) iJL{dz) . Observe that the kernel g must now satisfy 
/o /o '^^'^{^Tug(z))Ii{du)iJL{dz) < 00 

Assume now for simplicity that H{du) = H{u)du for some function 11 and that g is a non-negative 
monotone function with inverse g~^. In this case, we can obtain the LK form corresponding to II'- 



.-\lL(g) _ 



cxpf-/ / (l-e-^"9(^))n' 



{u)dufi{dz) 



exp 



il-e-'nKi^)dv 



where we have made the change of variables v — ug{z), z' ~ z and the measure Ilg{dv) is given by 

1 



n;(-) 



n' 



9iz')J gi^') 



j-fJ,{dz'). 



(32) 
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To demonstrate our method in this case, consider 



Il{9) 



-z/7-1 



L{dz 



00 /"OO 



ue-''^'^N{du,dz) 



where g{z) = e~^^^ with 77 > 0, // Lebesgue, and Il{du) = nu^^e^^^^du, which is the Levy measure 
corresponding to the Gamma distribution with shape k > and scale > ([3], Example 1.3.22). With the 
change of variables w = ue^^"^, ([32]) implies 



n'Ju) = / n' ue^/" e^/^'dz = ^ / w-'e-'^/'dw = ^T 0, ^ 



T^-n fr, " 



We can now compute the corresponding Laplace exponent for this case: 

^{X) = (1 - er^'')ii'(u)du = 77K / ^ — — ^r (0, -^ ) du 



u 

CXD pX pOO 

ut 



= T]K e 

'0 Jo 
r-x r 



(w-^e~'"/'>'\ dw dt du 



rjK 



rjK 



rjK 



rjK 



00 pW 



(w-ie-"*-"'/^)du dw 



dt 
dt 



'-Jo 
1 ' 
1 



-wjb 



./o 



-»(^+i/e)d^ ds 



rfs dwdt 



dt 



„.;'"■«"; "''di = „.L.(l + A9) 



^0 



The derivatives of c/) can also be computed exactly in this case. Using (|lip and ((2T|) . 
0(")(A) = (-l)"+i?7K /" M"-ie-^"r fo, ^"j du 

/■OO /"CXD 

-l)"+i77K / u"~ie-^" / (w-^e-"'/'')d«;dM 

Jo Ju 



io 

/•oo 

"'0 



^ ' A" 
^ A" 



^„-lg-An^^ 



dit; 



(n-l)!_^_,„,g(n-l)! ^" 



A" 



-1 -t«/6 



i! A'^ 



dw 



(1-e 



71 1 A^TT, /.C 



(33) 



(34) 



w 



^m-l^-^iX+l/e)^^ 



log(l + A0)-5^^ 
^^ m 1 



(A0)' 



'^ m (1 + xey 
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where the first integral in the second to last line above is computed as in (P^]) . Alternatively, ^^"•'(A) above 
can also be given simply in terms of the Gauss hypergeometric function 2^1 (a, /?, 7, x) ([9], Equation 6.455.1, 
page 657) 



^)(A) = (-1) 



rjK- 



2Fi(n,n,n + l,A6'). 



In Figure [7] we've plotted the PDF and CDF of (pij) for 6* = 1 and various values of the product p ~ tjk 
using our method. 
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Figure 7: Plots of the PDF and CDF for the Poisson stochastic integral /L(e */'') with various values of 
the parameter 77. 



6 Guide to software 

In this section, we will explain how to use the software written to implement the method discussed in this 
paper. Versions of this code exist in MATLAB and Mathematica, and are freely available by request form 
the authors. Each version will include a file containing examples to assist in using the code. 

To begin using the code, download the file NNINFDIV.zip and extract the directory NNINFDIV. This 
directory contains both the MATLAB and Mathematica programs in separate folders. We will now focus on 
these separately. 

6.1 MATLAB version 

To use the MATLAB version, launch MATLAB and add the directory NNINFDIV/MATLAB to MATLAB's 
default path by typing 

» path (path, 'mypath/NNINFDIV/MATLAB' ) 
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where 'mypath' is the path which leads to the directory NNINFDIV. You arc now ready to use the code 
provided in this package. 

The main function is called nninf div. This function takes in 5 arguments in the following order: 

X Scalar or vector of input values. Must be positive. 

DIST A cell array which specifies the distribution and parameters (see below) 

FUNC Type of function: 'pdf or 'cdf 

METHOD Extrapolation method: 'polynomial' or 'rational' 

TDL Target relative error tolerance 

The last three inputs FUNC, METHOD and TOL are optional, and take default values 'pdf , 'polynomial' and 
10^^. The input DIST is a cell array which contains the name of the desired distribution followed by the 
parameters. Possibilities for DIST are 

{ ' chi-squared' ,df , 1} Chi-squared distribution with df degrees of freedom. 

{ ' chi-squared' ,df , [cl, . . . ,cn]} Weighted sum of chi-squared distributions 

{'alpha stable' ,a,c} Alpha-stable distribution with a = a and scaling c. 

{ ' alpha stable ' , [al , . . . , an] , [cl , . . . , en] } Sum of n weighted alpha-stable distributions. 

{'uniform mix'} Uniform mix from Section 15.1.31 

{'ou poisson' ,eta} The integral /(e"''/'') from Scction r5.2.1l 

{'ou gamma' , eta, kappa} The integral /L(e~'*/'') from Section fS. 2. 21 

The scaling constants seen in the alpha-stable and chi-squared examples above compute the PDF/CDF of 
the scaled random variables cX for c > in the single alpha-stable and single chi-squared case. Likewise, in 
the weighted case, they return the PDF/CDF of ciXi + C2X2 + ■ ■ ■ c„X„ with q > 0, and Xi is chi-squared 
with df degrees of freedom or is alpha-stable with a =ai. 

Since the input for nninf div is long, it is often useful to define a function handle in order to call the 
function more easily. For example, consider the a-stable distribution with a = 2/3. We define the PDF of 
this distribution in the variable f by typing 

>f = @(x) nninfdiv(x,{ 'alpha stable' ,2/3,1}, 'pdf ', 'polynomial' ,le-6) ; 

The function f now computes the PDF of the alpha-stable distribution with a ~ 2/3 to within a relative 
error of 10^^ using the polynomial interpolation method. For example, you may now type 

> f (1) Computes the PDF at a; = 1 

> f ( [1 2 3] ) Computes the PDF at a; = 1, 2 and 3. 
» plot ( [ . 05 : . 05 : 2] , f ( [ . 05 : . 05 : 2] ) ) Plots the PDF on the interval (0,2] 

Remark: Obtaining relative errors less than 10^^ is sometimes difficult. If your error tolerance cannot be 
reached, the program will return the best estimate possible in double precision. If high precision is preferred 
over speed, the Mathematica version should be used. 

6.2 Mathematica version 

To use the Mathematica version, launch Mathematica open the file NNINFDIV. nb located in the directory 
NNINFDIV/Mathematica. Once this file is open, select all its contents by pressing alt-a on a PC or cmd-a on 
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a Mac. Then compile the code by pressing shift-return. You are now ready to use the code in a separate 
notebook. 

The main program is called NNInfDiv (capitalization matters). This program is called with 4 arguments: 

X Input value. Must be a positive scalar. 

DIST A list which specifies the distribution and parameters (see below) 

FUNC Type of function: "PDF" or "CDF" 

TOL Relative error tolerance 

The last two inputs FUNC and TOL are optional, taking default values "PDF" and 10"^'"' respectively. Possi- 
bilities for DIST include 

{"Chi-Squared" ,{cl , . . en}} Sum of weighted chi-squared with weights cl , . . . ,cn. 

{"Alpha Stable" ,a,c} Alpha-stable distribution with a = a and scaling c. 

{"Uniform Mix"} Uniform mix from Section fS. 1.31 

{"OU Poisson" ,eta} The integral /(e"^/'') from Section [?XT] 

{"OU Gamma" , eta, kappa} The integral Ihie'"^"^) from Section [021 

The scaling constants seen in the chi-squared and alpha-stable cases above refer to the random variables cX 
with c > in the alpha-stable case and ciXi -\- . . .CnXn in the chi-squared case, with Xi i.i.d chi-squared. 

To simplify the call to this function, one can make a user defined function. For example, to make a 
function F which computes the CDF of an a-stable distribution with a — 2/3, one can type 

F[xJ := NNlnfDiv[x,{ "Alpha Stable" ,2/3,1}, "CDF" ] 

The function F now computes the CDF of a-stable distribution with a — 2/3 to a relative precision of 10~^^. 
For example, one can now enter 

F [1] Computes the CDF at a: = 1 

Table [F[x] , {x, {1 ,2,3}}] Computes the CDF at a: = 1,2 and 3 

ListPlot [Table [{x,F[x]},{x, 0,2, .05}] , Joined -> True] Plots the CDF on the interval [0,2] 

Remark: Using NNInfDiv with Mathematica's Plot function is very slow, which is why we used the 
ListPlot function above. For faster plotting and function evaluation, the MATLAB version of the code 
should be used. 
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